Phenology



1. Load packages

require(ggplot2)
require(tidyverse)
require(dplyr)
require(nlme)
require(viridis)
require(data.table)
require(car)
require(lubridate)
require(stringr)
require(broom)


2. Load data

# define paths
data.dir   = "/Users/consti/Desktop/PhD/Publication_material/10_PhenogrowthEx3/Rdata/Experiment"
output.dir = "/Users/consti/Desktop/PhD/Publication_material/10_PhenogrowthEx3/Routput"

#load data
pheno      = data.frame(fread(paste(data.dir,"phenology.csv",sep="/")))         # Phenology data
biomass    = data.frame(fread(paste(data.dir,"biomass.csv",sep="/")))           # Biomass data
volume     = data.frame(fread(paste(data.dir,"volume.csv",sep="/")))            # Stem volume data

#get RMFs
biomass$RSR = biomass$Root/biomass$Shoot

#delete species epithet
DeleteEpithet <- function(df) {
    df$Species = gsub(' [A-z ]*', '' , df$Species)
    return(df)
}
dfnames    = list(pheno,biomass,volume)
dfs        = lapply(dfnames, DeleteEpithet)
names(dfs) = c("pheno","biomass","volume")

#Ordering of factors
Treatments      = c("control", "autumn", "autumn-spring", "spring") #order treatments
pheno$Treatment = factor(pheno$Treatment, levels=Treatments, ordered=T)
Species         = c("Quercus", "Fagus", "Lonicera") #order species
pheno$Species   = factor(pheno$Species, levels=Species, ordered=T)

#define plot themes
plotTheme = theme(legend.position   = "right",
                  legend.background = element_rect(fill=NA, size=0.5, linetype="solid"),
                  legend.text       = element_text(color="black"),
                  panel.grid.major  = element_blank(),
                  panel.grid.minor  = element_blank(),
                  panel.background  = element_rect(colour = NA, size=1, fill="grey97"),
                  axis.text.x       = element_text(angle = 45, hjust = 1),
                  axis.line.y       = element_line(color = "black"),
                  strip.background  = element_rect(fill=NA),
                  strip.text        = element_text(colour = 'black'))
plotTheme2 = theme(legend.position   = "right",
                  legend.background = element_rect(fill=NA, size=0.5, linetype="solid"),
                  legend.text       = element_text(color="black"),
                  panel.grid.major  = element_blank(),
                  panel.grid.minor  = element_blank(),
                  panel.background  = element_rect(colour = NA, size=1, fill="grey97"),
                  axis.line.y       = element_line(color = "black"),
                  strip.background  = element_rect(fill=NA),
                  strip.text        = element_text(colour = 'black'))


3. Prepare phenology data

#reshape phenology pheno to long format
pheno       = melt(dfs$pheno, id.vars = c("Species","ID","Treatment"))

#Convert dates to DOY
pheno$value = dmy(pheno$value)
pheno$DOY   = yday(pheno$value) # Jan 1 = day 1 

#create year and phenophase columns
pheno$Year  = str_sub(pheno$variable, str_length(pheno$variable)-3, str_length(pheno$variable))
pheno$Pheno = str_sub(pheno$variable, 1, str_length(pheno$variable)-5)
pheno       = dplyr::select(pheno, -c(variable, value))
pheno       = na.omit(pheno)

#from long (many rows) to short (multiple columns) format
pheno       = dcast(pheno, Species + ID + Treatment + Year ~ Pheno, value.var = "DOY")

#get growing season length information (leaf-out to senescence)
pheno$GSL   = pheno$Senescence - pheno$Leaf_out

#constrain late leafers
pheno$Leaf_out = ifelse(pheno$Leaf_out>140,140,pheno$Leaf_out)


4. Analysis


4.1. Treatment efffect on phenology

#################################
# Summarize annual phenology data
#################################


#get average phenology anomalies per individual
###############################################

#reshape phenology columns to long format
data = melt(pheno, id.vars = c("Species","ID","Treatment","Year"), 
             variable.name="phenology", value.name="day")

#get means per species and year
phenoMean = data %>% 
    group_by(Species,phenology,Year) %>% 
    summarize(Mean = mean(day, na.rm=TRUE))

#merge data
data <-merge(data, phenoMean, all.x=T, by=c("Species","phenology","Year"))
rm(phenoMean)

#get change in phenology relative to mean
data = data %>% 
  mutate(phenologyAnomaly = ifelse(phenology=="Leaf_out", Mean-day, day-Mean)) %>% #add column
  dplyr::select(!c(Mean)) #delete columns

#from long (many rows) to short (multiple columns) format
data3 = dcast(data, Species + ID + Treatment + Year ~ phenology, value.var = c("phenologyAnomaly"))

#summarize years
VariablesVector = c("Leaf_out","Senescence","GSL")
data3 = data.frame(data3 %>%
    filter(!Year %in% c("2017")) %>% #delete rows based on condition
    group_by(Species,ID,Treatment) %>% 
    summarize_at(VariablesVector, mean, na.rm = TRUE))


# Add biomass information
#########################

#merge phenology and biomass tables
data3 = merge(data3, dfs$biomass, all.x=T, by=c("Species","ID","Treatment"))

#reshape phenology columns to long format
data3 = melt(data3, id.vars = c("Species","ID","Treatment","Total","Root","Shoot","RSR"), 
             variable.name="phenology", value.name="phenologyAnomaly")


#get average phenology dates per individual
###########################################

#from long (many rows) to short (multiple columns) format
data1 = dcast(data, Species + ID + Treatment + Year ~ phenology, value.var = c("day"))

#summarize years
VariablesVector = c("Leaf_out","Senescence","GSL")
data1 = data.frame(data1 %>%
    filter(!Year %in% c("2017")) %>% #delete rows based on condition
    group_by(Species,ID,Treatment) %>% 
    summarize_at(VariablesVector, mean, na.rm = TRUE))

#merge phenology and biomass tables
data1 = merge(data1, dfs$biomass, all.x=T, by=c("Species","ID","Treatment"))

#reshape phenology columns to long format
data1 = melt(data1, id.vars = c("Species","ID","Treatment","Total","Root","Shoot","RSR"), 
             variable.name="phenology", value.name="day")


#########################################################################################
#########################################################################################


###########################
# Multivariate linear model
###########################


#Create spring and autumn treatments
####################################

data4        = data3
data4$spring = ifelse(data4$Treatment %in% c("control","autumn"), "no","warming")
data4$autumn = ifelse(data4$Treatment %in% c("control","spring"), "no","warming")

#Multivariate linear model
##########################

resultsLM = data4 %>% 
    group_by(Species,phenology) %>% 
    do({model = lm(phenologyAnomaly ~ spring*autumn, data=.)  # create your model
    data.frame(tidy(model))}) %>%           # get coefficient info
    filter(!term %in% c("(Intercept)")) %>% # delete rows
    dplyr::select(Species, phenology, term, p.value, estimate, std.error) %>% #delete columns
    mutate_if(is.numeric, round, 2) %>%   
    mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) #add significance
#Ordering of factors
term           = c("springwarming", "autumnwarming", "springwarming:autumnwarming") #order treatments
resultsLM$term = factor(resultsLM$term, levels=term, ordered=T)
Species         = c("Quercus", "Fagus", "Lonicera") #order species
resultsLM$Species   = factor(resultsLM$Species, levels=Species, ordered=T)
#order statistics table to match dataframe
resultsLM = resultsLM[with(resultsLM, order(resultsLM$Species, resultsLM$phenology)), ]
data.frame(resultsLM[,-c(7)])
##     Species  phenology                        term p.value estimate std.error
## 1   Quercus   Leaf_out               springwarming    0.00    16.27      1.16
## 2   Quercus   Leaf_out               autumnwarming    0.06    -2.17      1.09
## 3   Quercus   Leaf_out springwarming:autumnwarming    0.79     0.43      1.57
## 4   Quercus Senescence               springwarming    0.02   -10.77      4.56
## 5   Quercus Senescence               autumnwarming    0.00    13.67      4.31
## 6   Quercus Senescence springwarming:autumnwarming    0.39    -5.43      6.19
## 7   Quercus        GSL               springwarming    0.27     5.50      4.88
## 8   Quercus        GSL               autumnwarming    0.02    11.50      4.61
## 9   Quercus        GSL springwarming:autumnwarming    0.46    -5.00      6.63
## 10    Fagus   Leaf_out               springwarming    0.00    13.72      1.66
## 11    Fagus   Leaf_out               autumnwarming    0.00    -6.73      1.66
## 12    Fagus   Leaf_out springwarming:autumnwarming    0.71     0.89      2.37
## 13    Fagus Senescence               springwarming    0.12    -3.28      2.05
## 14    Fagus Senescence               autumnwarming    0.00    23.11      2.05
## 15    Fagus Senescence springwarming:autumnwarming    0.74    -1.00      2.94
## 16    Fagus        GSL               springwarming    0.00    10.44      2.21
## 17    Fagus        GSL               autumnwarming    0.00    16.38      2.21
## 18    Fagus        GSL springwarming:autumnwarming    0.97    -0.11      3.16
## 19 Lonicera   Leaf_out               springwarming    0.00    21.88      1.95
## 20 Lonicera   Leaf_out               autumnwarming    0.09    -3.54      2.06
## 21 Lonicera   Leaf_out springwarming:autumnwarming    0.00    -8.84      2.83
## 22 Lonicera Senescence               springwarming    0.04    -4.75      2.19
## 23 Lonicera Senescence               autumnwarming    0.00     7.69      2.32
## 24 Lonicera Senescence springwarming:autumnwarming    0.08     5.67      3.19
## 25 Lonicera        GSL               springwarming    0.00    17.13      2.77
## 26 Lonicera        GSL               autumnwarming    0.17     4.15      2.93
## 27 Lonicera        GSL springwarming:autumnwarming    0.44    -3.17      4.04
#Plot
#####

#create panel labels
dat_text       = data.frame(
  label        = c("a","b","c","d","e","f","g","h","i"),
  phenology   = rep(c("Leaf_out","Senescence","GSL"),3),
  Species      = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
  term         = "springwarming")
#Plot
LMplot = ggplot(data = resultsLM, mapping = aes(x = term, y = estimate, 
                                                fill=phenology, alpha=term)) +
  geom_bar(position=position_dodge(), stat="identity", color="black") +
    geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+
  geom_hline(yintercept=0) +
  geom_text(data = dat_text, mapping = aes(x = Inf, y = Inf, 
                                           hjust = 1, vjust = 1,
                                           label = paste("(",label,")",sep=""),
                                           fontface = "bold"))+
  stat_summary(geom = 'text', label = resultsLM$sig, 
               fun.y = mean, vjust = 2) +
  scale_alpha_discrete(range = c(1, 0.7)) +
  xlab("Phenophase") +
  ylab("Effect size (days)") +
  scale_fill_manual(values = c("green3", "orange", "blue")) +
  facet_grid(Species~phenology) + 
  plotTheme
LMplot

######################################################
# Check for deviation from control (one-sample T-test)
######################################################

#get means of control group
data2 = data1 %>% 
    filter(Treatment %in% c("control")) %>%
    group_by(Species,phenology) %>% 
    summarize(Mean = mean(day, na.rm=TRUE))
#merge data
data2 <-merge(data1[!data1$Treatment=="control",], data2, all.x=T, by=c("Species","phenology"))
#get change in phenology relative to control
data2$phenologyChange = data2$day - data2$Mean

#Ordering of factors
Phenology     = c("Leaf_out", "Senescence", "GSL") #order treatments
data2$phenology = factor(data2$phenology, levels=Phenology, ordered=T)
Treatments     = c("spring", "autumn", "autumn-spring") #order treatments
data2$Treatment = factor(data2$Treatment, levels=Treatments, ordered=T)
Species         = c("Quercus", "Fagus", "Lonicera") #order species
data2$Species   = factor(data2$Species, levels=Species, ordered=T)

#T-test
resultsTT = data2 %>% 
    group_by(Species,phenology,Treatment) %>% 
    do({model = t.test(.$phenologyChange) # create your model
    data.frame(tidy(model),
               tidy(shapiro.test(x = .$phenologyChange))[1,2])}) %>%           
    dplyr::select(Species, phenology, Treatment, p.value, estimate, p.value.1) %>% #delete columns
    rename(shapiroTest=p.value.1) %>%      
    mutate(significance = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) %>% #add significance
    mutate_if(is.numeric, round, 2) 
#order statistics table to match dataframe
resultsTT = resultsTT[with(resultsTT, order(resultsTT$Species, resultsTT$phenology)), ]
data.frame(resultsTT[,-c(7)])
##     Species  phenology     Treatment p.value estimate shapiroTest
## 1   Quercus   Leaf_out        spring    0.00   -16.27        0.41
## 2   Quercus   Leaf_out        autumn    0.01     2.17        0.65
## 3   Quercus   Leaf_out autumn-spring    0.00   -14.53        0.24
## 4   Quercus Senescence        spring    0.00   -10.77        0.36
## 5   Quercus Senescence        autumn    0.00    13.67        0.73
## 6   Quercus Senescence autumn-spring    0.52    -2.53        0.37
## 7   Quercus        GSL        spring    0.03     5.50        0.33
## 8   Quercus        GSL        autumn    0.01    11.50        0.69
## 9   Quercus        GSL autumn-spring    0.01    12.00        0.38
## 10    Fagus   Leaf_out        spring    0.00   -13.72        0.41
## 11    Fagus   Leaf_out        autumn    0.00     6.73        0.17
## 12    Fagus   Leaf_out autumn-spring    0.00    -7.88        0.77
## 13    Fagus Senescence        spring    0.08    -3.28        0.49
## 14    Fagus Senescence        autumn    0.00    23.11        0.02
## 15    Fagus Senescence autumn-spring    0.00    18.83        0.22
## 16    Fagus        GSL        spring    0.00    10.44        0.60
## 17    Fagus        GSL        autumn    0.00    16.38        0.31
## 18    Fagus        GSL autumn-spring    0.00    26.72        0.42
## 19 Lonicera   Leaf_out        spring    0.00   -21.88        0.08
## 20 Lonicera   Leaf_out        autumn    0.01     3.54        0.36
## 21 Lonicera   Leaf_out autumn-spring    0.00    -9.50        0.98
## 22 Lonicera Senescence        spring    0.02    -4.75        0.48
## 23 Lonicera Senescence        autumn    0.00     7.69        0.02
## 24 Lonicera Senescence autumn-spring    0.00     8.61        0.01
## 25 Lonicera        GSL        spring    0.00    17.13        0.44
## 26 Lonicera        GSL        autumn    0.00     4.15        0.36
## 27 Lonicera        GSL autumn-spring    0.00    18.11        0.98
# Plot
######

#create panel labels
dat_text <- data.frame(
  label       = c("a","b","c","d","e","f","g","h","i"),
  phenology   = rep(c("Leaf_out","Senescence","GSL"),3),
  Species     = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
  Treatment   = "spring")
#Plot
PhenologyChangePlot = ggplot(data = data2, mapping = aes(x = Treatment, y = phenologyChange, 
                                                         fill=phenology, alpha=Treatment)) +
  stat_summary(fun.y = mean, 
                 geom = "bar", color="black") + 
  stat_summary(fun.data = mean_se,  
                 geom = "errorbar", width=0.2) + 
  geom_hline(yintercept=0) +
  geom_text(data = dat_text, mapping = aes(x = -Inf, y = Inf, 
                                           hjust = -0.2, vjust = 1.5,
                                           label = paste("(",label,")",sep=""),
                                           fontface = "bold"))+
  stat_summary(geom = 'text', label = resultsTT$significance, 
               fun.y = mean, vjust = 2) +
  scale_alpha_discrete(range = c(1, 0.7)) +
  xlab("Treatment") +
  ylab("Phenological change (days)") +
  scale_fill_manual(values = c("green3", "orange", "blue")) +
  facet_grid(Species~phenology) + 
  plotTheme
PhenologyChangePlot

##########
#Save PDFs
##########

pdf(paste(output.dir,"PhenologyLinearTreatmentModel.pdf",sep="/"), width=7, height=7, useDingbats=FALSE)
LMplot
dev.off()
## quartz_off_screen 
##                 2
pdf(paste(output.dir,"PhenologyChange.pdf",sep="/"), width=7, height=7, useDingbats=FALSE)
PhenologyChangePlot
dev.off()
## quartz_off_screen 
##                 2


4.2. Phenology effect on total biomass

##########################
# Univariate linear models
##########################

#reshape biomass columns to long format
data3 = melt(data3, id.vars = c("Species","ID","Treatment","phenology","phenologyAnomaly"), 
             variable.name="organ", value.name="biomass")

#Transform biomass to percentage anomaly
########################################

#get means per species
data5 = data3 %>% 
    group_by(Species,organ) %>% 
    summarize(Mean = mean(biomass, na.rm=TRUE))
#merge data
data3 <-merge(data3, data5, all.x=T, by=c("Species","organ"))
#get precentage change in biomass
data3$relativeBiomass = (data3$biomass / data3$Mean -1) * 100

#order species
Species           = c("Quercus", "Fagus", "Lonicera") 
data3$Species = factor(data3$Species, levels=Species, ordered=T)

#Extract linear model info
resultsLM = data3 %>% 
    group_by(Species,phenology,organ) %>% 
    do({model = lm(relativeBiomass ~ phenologyAnomaly, data=.)  # create model
    data.frame(tidy(model),                # get coefficient info
               glance(model),
               tidy(shapiro.test(x = residuals(object = model)))[1,2]
               )}) %>%                     # get model info
    filter(term != "(Intercept)") %>%      # delete intercept info
    dplyr::select(Species, phenology, organ, estimate, std.error, p.value, r.squared, p.value.2) %>% # delete columns
    rename(shapiro=p.value.2) %>%          #rename columns
    mutate(sig = ifelse(p.value<0.05, "**", ifelse(p.value<=0.1, "*", ""))) %>% # add column
    mutate_if(is.numeric, round, 2)

#Order phenophases
Organs            = c("Total","Shoot","Root","RSR") #order treatments
resultsLM$organ   = factor(resultsLM$organ, levels=Organs, ordered=T)
#order statistics table to match dataframe
resultsLM = resultsLM[with(resultsLM, order(resultsLM$Species, resultsLM$organ, resultsLM$phenology)), ]
data.frame(resultsLM)
##     Species  phenology organ estimate std.error p.value r.squared shapiro sig
## 1   Quercus   Leaf_out Total     2.46      0.77    0.00      0.22    0.39  **
## 2   Quercus Senescence Total    -2.12      0.49    0.00      0.35    0.95  **
## 3   Quercus        GSL Total    -1.30      0.67    0.06      0.10    0.35   *
## 4   Quercus   Leaf_out Shoot     1.79      0.74    0.02      0.14    0.23  **
## 5   Quercus Senescence Shoot    -1.98      0.44    0.00      0.37    0.61  **
## 6   Quercus        GSL Shoot    -1.53      0.59    0.01      0.16    0.22  **
## 7   Quercus   Leaf_out  Root     3.03      0.89    0.00      0.25    0.62  **
## 8   Quercus Senescence  Root    -2.25      0.60    0.00      0.29    0.75  **
## 9   Quercus        GSL  Root    -1.10      0.81    0.18      0.05    0.02    
## 10  Quercus   Leaf_out   RSR     0.98      0.54    0.08      0.09    0.83   *
## 11  Quercus Senescence   RSR    -0.18      0.39    0.64      0.01    0.16    
## 12  Quercus        GSL   RSR     0.38      0.45    0.40      0.02    0.36    
## 13    Fagus   Leaf_out Total     0.78      0.51    0.14      0.06    0.28    
## 14    Fagus Senescence Total    -0.16      0.36    0.66      0.01    0.60    
## 15    Fagus        GSL Total     0.27      0.41    0.52      0.01    0.44    
## 16    Fagus   Leaf_out Shoot     0.42      0.52    0.43      0.02    0.31    
## 17    Fagus Senescence Shoot     0.09      0.36    0.80      0.00    0.34    
## 18    Fagus        GSL Shoot     0.37      0.40    0.36      0.02    0.50    
## 19    Fagus   Leaf_out  Root     1.24      0.60    0.05      0.11    0.18  **
## 20    Fagus Senescence  Root    -0.49      0.43    0.27      0.03    0.57    
## 21    Fagus        GSL  Root     0.13      0.50    0.80      0.00    0.07    
## 22    Fagus   Leaf_out   RSR     0.81      0.42    0.06      0.10    0.28   *
## 23    Fagus Senescence   RSR    -0.50      0.29    0.10      0.08    0.82    
## 24    Fagus        GSL   RSR    -0.14      0.34    0.68      0.00    0.94    
## 25 Lonicera   Leaf_out Total    -0.16      0.20    0.45      0.02    0.39    
## 26 Lonicera Senescence Total     0.37      0.30    0.22      0.04    0.44    
## 27 Lonicera        GSL Total     0.01      0.22    0.96      0.00    0.44    
## 28 Lonicera   Leaf_out Shoot    -0.13      0.22    0.58      0.01    0.75    
## 29 Lonicera Senescence Shoot     0.07      0.34    0.84      0.00    0.90    
## 30 Lonicera        GSL Shoot    -0.12      0.25    0.64      0.01    0.79    
## 31 Lonicera   Leaf_out  Root    -0.21      0.30    0.50      0.01    0.50    
## 32 Lonicera Senescence  Root     0.91      0.43    0.04      0.12    0.09  **
## 33 Lonicera        GSL  Root     0.24      0.33    0.48      0.02    0.18    
## 34 Lonicera   Leaf_out   RSR    -0.09      0.33    0.79      0.00    0.27    
## 35 Lonicera Senescence   RSR     0.87      0.47    0.07      0.09    0.43   *
## 36 Lonicera        GSL   RSR     0.37      0.36    0.32      0.03    0.26
#Visually inspect model assumptions 
data.assumptions = data3
data.assumptions$category = paste(data.assumptions$Species,
                                  data.assumptions$organ,
                                  data.assumptions$phenology, sep="_") #create identifier column
category.list = as.factor(unique(data.assumptions$category))       #create category vector
par(mfrow=c(2,2))                                       #set plot layout
for (category in category.list){                        #loop over categories
  tab.subset=data.assumptions[data.assumptions$category==category, ]          #create table subset
  plot(lm(relativeBiomass ~ phenologyAnomaly, data=tab.subset), main=category)
}

######
#Plots
######

#Effect sizes
#############

#create panel labels
dat_text      = data.frame(
  label       = c("a","b","c","d","e","f","g","h","i","j","k","l"),
  organ       = rep(c("Total","Shoot","Root","RSR"),3),
  Species     = c(rep("Quercus",4), rep("Fagus",4), rep("Lonicera",4)),
  phenology   = "Leaf_out")

#Plot
LMplot = ggplot(data = resultsLM, mapping = aes(x = phenology, y = estimate, 
                                                fill=organ, alpha=phenology)) +
  geom_bar(position=position_dodge(), stat="identity", color="black") +
    geom_errorbar(aes(ymin=estimate-std.error, ymax=estimate+std.error),
                  width=.2,                    # Width of the error bars
                  position=position_dodge(.9))+
  geom_hline(yintercept=0) +
  geom_text(data = dat_text, mapping = aes(x = -Inf, y = -Inf, 
                                           hjust = -0.1, vjust = -1,
                                           label = paste("(",label,")",sep=""),
                                           fontface = "bold"))+
  stat_summary(geom = 'text', label = resultsLM$sig, 
               fun.y = mean, vjust = 2) +
  scale_color_viridis(option="viridis",discrete=TRUE) +
  scale_alpha_discrete(range = c(1, 0.7)) +
  xlab("Phenophase") +
  ylab("Effect size (%/day)") +
  scale_fill_viridis(option="E",discrete=TRUE) +
  facet_grid(Species~organ) + 
  plotTheme
LMplot

#Biomass
########

#order
Organs      = c("Total","Shoot","Root","RSR") #order treatments
data3$organ = factor(data3$organ, levels=Organs, ordered=T)
#create panel labels
dat_text       = data.frame(
  label        = c("a","b","c","d","e","f","g","h","i"),
  phenology    = rep(c("Leaf_out","Senescence","GSL"),3),
  Species      = c(rep("Quercus",3), rep("Fagus",3), rep("Lonicera",3)),
  organ        = "Root")
BiomassPlot = ggplot(data3[data3$organ %in% c("Total","Shoot","Root"),], 
                     aes(x=phenologyAnomaly, y=relativeBiomass, colour=organ)) + 
    #geom_point(size=0.2) +  
    geom_hline(yintercept=0) +
    geom_vline(xintercept=0) +
    geom_smooth(method=lm, fullrange = T, alpha = 0.2) +
    labs(x = "Phenology anomaly (days)", y = "Biomass anomaly (%)") +
    scale_color_manual(values=c("orange","yellow2","red3"))+
    geom_text(data = dat_text, mapping = aes(x = -Inf, y = Inf, 
                                           hjust = -0.2, vjust = 1.3,
                                           label = paste("(",label,")",sep=""),
                                           fontface = "bold"), color="black")+
    geom_text(data    = resultsLM[resultsLM$organ=="Total",],
              mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 2.5, 
                            label = paste(round(estimate,1), " %/day", sig, sep="")), 
              size=3.5, color="orange")+
    geom_text(data    = resultsLM[resultsLM$organ=="Shoot",],
              mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 4.5, 
                            label = paste(round(estimate,1), " %/day", sig, sep="")), 
              size=3.5, color="yellow2")+
    geom_text(data    = resultsLM[resultsLM$organ=="Root",],
              mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 6.5, 
                            label = paste(round(estimate,1), " %/day", sig, sep="")), 
              size=3.5, color="red3")+
    facet_grid(Species~phenology, scales = "free") +
    plotTheme2
BiomassPlot

#RSR
RSRplot = ggplot(data3[data3$organ %in% c("RSR"),], aes(x=phenologyAnomaly, y=biomass)) + 
    geom_point(size=0.2) +  
    geom_smooth(method=lm, fullrange = F) +
    labs(x = "Day", y = "Root:shoot ratio") +
    scale_color_viridis(option="viridis",discrete=TRUE) +
    geom_text(data    = resultsLM[resultsLM$organ=="RSR",],
              mapping = aes(x = Inf, y = Inf, hjust = 1.5, vjust = 2.5, 
                            label = paste("Slope = ", round(estimate,3), " day-1", sig, sep="")), 
              size=3.5, color="black")+
    facet_grid(Species~phenology, scales = "free") +
    plotTheme2 
RSRplot

#########
#Save PDF
#########

pdf(paste(output.dir,"PhenophasesEffectSizes.pdf",sep="/"), width=8, height=7, useDingbats=FALSE)
LMplot
dev.off()
## quartz_off_screen 
##                 2
pdf(paste(output.dir,"AveragePhenology.pdf",sep="/"), width=8, height=7, useDingbats=FALSE)
BiomassPlot
RSRplot
dev.off()
## quartz_off_screen 
##                 2